Miopia ou controle? Análise da eficiência dos instrumentos regulatórios urbanos para determinação de densidade habitacional

library(tidyverse)

library(sf)
library(exactextractr)

library(patchwork)
library(ggtext)
library(units)

library(modelsummary)
library(gt)


set.seed(123)

distrito <- read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% st_set_crs("epsg:31983")

Microdados do Censo de 2022

É importante destacar que os resultados disponíveis são parciais do censo, com dados atualizados até o momento da divulgação. No entanto, mesmo parciais, esses dados podem ser extremamente úteis para suas análises.

Um ponto crucial a considerar é que o nível mínimo de observação georreferenciada nos microdados do censo são os setores censitários. Os setores censitários são unidades geográficas definidas pelo Instituto Brasileiro de Geografia e Estatística (IBGE) para coleta e tabulação de dados censitários. Eles são delimitados de forma a garantir uma cobertura completa e homogênea de todo o território nacional, facilitando a análise comparativa entre diferentes áreas geográficas. A delimitação geográfica do setor também considera questões logísticas para o entrevistador conseguir entrevistar a todos.

# Read dados do censo 2022
censo <- read_sf("dados/censo/SP_Malha_Preliminar_2022.shp") %>% 
  filter(CD_MUN == "3550308") %>% 
  st_transform("epsg:31983") %>% # Sistema de coordenadas do geosampa
  select(id_setor_censitario = CD_SETOR, v0001:v0007) %>% 
  mutate(area_setor = st_area(geometry) %>% as.numeric())

censo %>% 
  st_drop_geometry() %>% 
  modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
v0001	1225	0	415.0	230.9	0.0	397.0	2221.0
v0002	592	0	181.1	94.3	0.0	177.0	1372.0
v0003	588	0	180.9	94.3	0.0	177.0	1371.0
v0004	26	0	0.2	1.6	0.0	0.0	170.0
v0005	14492	0	2.6	0.6	0.0	2.7	10.3
v0006	7559	0	13.6	12.1	0.0	10.8	100.0
v0007	511	0	156.4	81.9	0.0	152.0	950.0
area_setor	27592	0	55126.5	360317.6	196.3	23142.2	22848344.5

gg <- ggplot() +
  geom_sf(data = censo %>%
            mutate(densidade = v0001/area_setor,
                   decil = ntile(densidade, 10)),
          aes(fill = factor(decil)), color = NA) +
  scale_fill_viridis_d() +
  labs(fill = "Decil de densidade") +
  theme_void() +
  theme(legend.position = c(.8, .3))

ggsave("tex/imagens/mapa.pdf", gg, width = 7, height = 12)
gg

censo %>% 
  st_drop_geometry() %>% 
  summarize("Total de Pessoas" = sum(v0001),
            "Total de Domicílios" = sum(v0002),
            "Total de Domicílios Particulares Ocupados" = sum(v0007)) %>% 
  pivot_longer(everything(), names_to = "Variável", values_to = "Valor") %>% 
  gt()

Variável	Valor
Total de Pessoas	11451999
Total de Domicílios	4996529
Total de Domicílios Particulares Ocupados	4316336

gg <- censo %>% 
  mutate(distancia_centro = st_distance(geometry, 
                                        read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% 
                                          st_set_crs("epsg:31983") %>% 
                                          filter(ds_nome == "SE")) %>% 
           as.numeric() %>% 
           cut(breaks = 10^3*(0:40), labels = FALSE)) %>% 
  st_drop_geometry() %>%
  group_by(distancia_centro) %>% 
  summarize(densidade = sum(v0001) / sum(as.numeric(area_setor) / 10^6)) %>% 
  mutate(teorico = 2*10^4*exp(-distancia_centro/10)) %>% 
  ggplot(aes(x = distancia_centro)) +
  geom_col(aes(y = densidade), fill = "#2BAA92FF") +
  geom_line(aes(y = teorico), lwd = .8, linetype = "dashed", color = "#D32934FF") +
  scale_y_continuous(labels = scales::comma_format(big.mark = ".")) +
  scale_x_continuous(labels = scales::comma_format(big.mark = ".")) +
  labs(x = "Distância ao centro em km", 
       y = expression("Densidade populacional por" ~ km^2)) +
  theme_classic()

ggsave("tex/imagens/densidade_distcentro.pdf", gg, width = 5, height = 3)
gg

Base de Dados do IPTU

A base de dados do IPTU (Imposto Predial e Territorial Urbano) é uma fonte abrangente de informações sobre imóveis urbanos dentro do município. Essa base é considerada completa, pois abrange todos os imóveis sujeitos à tributação do IPTU, representando uma fonte confiável e abrangente de informações sobre a propriedade urbana. Propriedades que não foram construídas dentro da legalidade não constam nessa base.

O número do contribuinte, utilizado para identificar exclusivamente cada imóvel na base de dados do IPTU, é representado diretamente pelo SQL, sendo essencial para consultas e manipulação dos dados relacionados ao imposto predial e territorial urbano. Esse formato permite a integração e cruzamento com outras bases de dados, como a base de lotes, que está georreferenciada, fornecendo informações espaciais adicionais que não estão disponíveis na base do IPTU.

Quando um lote é um condomínio, ou seja, quando não se classifica de acordo com condominio == "00-0", é necessário substituir os múltiplos números de SQL pelo número do condomínio. Isso ocorre porque cada unidade dentro do condomínio pode ser tratada como uma entidade separada para fins tributários, mas para a base de lotes, estão todos juntos.

IPTU <- read.csv("dados/IPTU/IPTU_2024.csv", sep=";", encoding = "latin1") %>% 
    as_tibble() %>% 
    select(sql = "NUMERO.DO.CONTRIBUINTE", 
           condominio = "NUMERO.DO.CONDOMINIO",
           area_terreno = "AREA.DO.TERRENO",
           area_construida = "AREA.CONSTRUIDA",
           area_ocupada = "AREA.OCUPADA",
           pavimentos = "QUANTIDADE.DE.PAVIMENTOS",
           ano_construcao = "ANO.DA.CONSTRUCAO.CORRIGIDO",
           tipo = "TIPO.DE.PADRAO.DA.CONSTRUCAO") %>% 
    
    # Separação do número de contribuinte (SQL) em setor quadra e lote
    mutate(setor =  str_sub(sql, 1, 3),
           quadra = str_sub(sql, 4, 6),
           
           # Quando o lote é um condomínio, haverá vários SQLs no mesmo lote. CD = Condomínio
           lote = str_sub(sql, 7, 10) %>% 
             ifelse(condominio == "00-0", ., paste("CD", str_sub(condominio, 1, 2), sep = "")),
           
           # Tipo de uso
           residencial = str_detect(tipo, "Residencial"))

IPTU %>% 
  modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
area_terreno	9304	0	3546.7	15146.5	1.0	800.0	4307493.0
area_construida	6718	0	157.1	1068.9	0.0	102.0	950328.0
area_ocupada	5159	0	1276.6	3170.0	0.0	378.0	320000.0
pavimentos	52	0	9.4	8.7	0.0	5.0	98.0
ano_construcao	127	0	1948.3	283.9	0.0	1987.0	2023.0

Todos os empreendimentos não públicos regulares (de acordo com a lei) constam nesta base de dados. A seguir está uma análise de quais são estes tipos. O tipo de interesse para esta análise é apenas os residenciais.

gg.right <- IPTU %>% 
  mutate(uso = case_when(str_detect(tipo, "Residencial") ~ "Residencial",
                         str_detect(tipo, "Comercial") ~ "Comercial",
                         str_detect(tipo, "Oficina") ~ "Industria",
                         str_detect(tipo, "TERRENO") ~ "Terreno",
                         str_detect(tipo, "Clube") ~ "Entretenimento",
                         .default = "Outros") %>% as.factor(),
         padrao = case_when(str_detect(tipo, "A$") ~ "A",
                            str_detect(tipo, "B$") ~ "B",
                            str_detect(tipo, "C$") ~ "C",
                            str_detect(tipo, "D$") ~ "D",
                            str_detect(tipo, "E$") ~ "E",
                            .default = "NA"),
         nome = paste("(", padrao, ") ", uso, sep = "")) %>% 
  group_by(uso, padrao, nome) %>% 
  summarize(area = sum(area_construida)) %>% 
  ungroup() %>% 
  mutate(percentual = area * 100/ sum(area, na.rm = TRUE),
         texto = case_when(percentual < 1 ~ "<1%",
                           percentual < 5 ~ paste(round(percentual,1), "%", sep = ""),
                           .default = paste(round(percentual, 2), "%", sep = "")),
         color = ifelse(uso == "Outros", "white", "black")) %>% 
  ggplot() +
  treemapify::geom_treemap(aes(fill = uso, area = area), size = 2, color = NA) +
  treemapify::geom_treemap(aes(alpha = padrao, area = area), fill = "black", color = NA) +
  treemapify::geom_treemap_text(aes(area = area, label = nome, color = color), alpha = .4, grow = FALSE, size = 8) +
  treemapify::geom_treemap_text(aes(area = area, label = texto, color = color), 
                                alpha = .4, place = "middle", size = 10) +
  scale_fill_manual(values = c("Residencial" = "#FAE48B", 
                               "Comercial" = "#A6EEE6",
                               "Industria" = "#84BD00",
                               "Entretenimento" = "#FB6467",
                               "Outros" = "#1A5354")) +
  scale_alpha_manual(values = c("A" = 0, "B" = .05, "C" = .1, "D" = .15, "E" = .2, "NA" = 0)) +
  scale_color_manual(values = c("black" = "black", "white" = "white")) +
  labs(fill = "Tipo de uso", alpha = "Padrão de uso") +
  theme(aspect.ratio=1)

gg.left <- IPTU %>% 
  mutate(uso = case_when(str_detect(tipo, "Residencial") ~ "Residencial",
                         str_detect(tipo, "Comercial") ~ "Comercial",
                         str_detect(tipo, "Oficina") ~ "Industria",
                         str_detect(tipo, "TERRENO") ~ "Terreno",
                         str_detect(tipo, "Clube") ~ "Entretenimento",
                         .default = "Outros") %>% as.factor()) %>% 
  group_by(uso) %>% 
  summarize(area = sum(area_construida, na.rm = TRUE)) %>% 
  ungroup() %>% 
  arrange(desc(uso)) %>% 
  mutate(percentual = area * 100 / sum(area, na.rm = TRUE),
         texto = ifelse(percentual > 5, paste(round(percentual, 1), "%", sep = ""), ""),
         ytext = (cumsum(area) + lag(cumsum(area)))/ 2) %>% 
  ggplot() +
  geom_col(aes(y = area, fill = uso, x = "")) +
  geom_text(aes(y = ytext, fill = uso, x = "", label = texto), alpha = .5) +
  scale_fill_manual(values = c("Residencial" = "#FAE48B", 
                               "Comercial" = "#A6EEE6",
                               "Industria" = "#84BD00",
                               "Entretenimento" = "#FB6467",
                               "Outros" = "#1A5354")) +
  theme_void() +
  theme(legend.position = "none")

ggsave("tex/imagens/tree_area_construida.pdf", 
       (gg.left | gg.right) + 
         plot_layout(widths = c(1,6)), 
       width = 6.75, height = 5)

(gg.left | gg.right)

Quando há um condomínio com múltiplos contribuintes de IPTU, ele deve ser agregado a nível lote, para que possa ser cruzado com a base de lotes, possibilitando que seja georreferenciada. Todos os contribuintes dentro de um condomínio compartilham das mesmas características de área de terreno, área ocupada, ano de construção e pavimentos, então a mediana funciona para agregar estes dados, mas o máximo, mínimo, média ou pegar o primeiro valor também funcionaria.

IPTU <- IPTU %>% 
  # Agregar por SQL
  group_by(setor, quadra, lote) %>% 
  summarize(unidades = n(),
            area_terreno = median(area_terreno), 
            area_construida = sum(area_construida), 
            area_ocupada = median(area_ocupada),
            pavimentos = median(pavimentos),
            ano_construcao = median(ano_construcao),
            residencial = median(residencial)) %>% 
  ungroup() %>% 
  mutate(CA = area_construida / area_terreno,
         cota_parte = area_terreno / unidades,
         verticalizacao = area_construida / area_ocupada) %>% 
  
  filter(area_ocupada > 0)

IPTU %>% 
  filter(residencial == 1) %>% 
  modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
unidades	545	0	2.4	16.4	1.0	1.0	2861.0
area_terreno	5736	0	248.8	1917.0	6.0	160.0	1442400.0
area_construida	12097	0	306.3	1954.9	2.0	120.0	400000.0
area_ocupada	3403	0	118.4	407.8	1.0	90.0	200000.0
pavimentos	51	0	1.8	1.8	1.0	2.0	52.0
ano_construcao	166	0	1981.8	14.6	1800.0	1981.0	2023.0
residencial	1	0	1.0	0.0	1.0	1.0	1.0
CA	100483	0	0.9	0.8	0.0	0.8	36.1
cota_parte	17359	0	208.3	1624.9	1.5	155.0	1442400.0
verticalizacao	53664	0	1.7	1.7	0.0	1.4	180.0

IPTU %>% 
  mutate(Tipo = ifelse(unidades == 1, "Unidade", "Condomínio")) %>% 
  group_by(Tipo) %>% 
  summarize("Lotes" = n(),
            "Unidades" = sum(unidades)) %>% 
  gt()

Tipo	Lotes	Unidades
Condomínio	33116	1781971
Unidade	1255235	1255235

gg <- IPTU %>% 
  filter(residencial == 1) %>% 
  filter(cota_parte < 600, CA < 5, verticalizacao < 6) %>% 
  select("Cota Parte" = cota_parte, "Coeficiente de Aproveitamento (CA)" = CA, "Verticalização" = verticalizacao) %>% 
  pivot_longer(everything()) %>% 
  ggplot() +
  geom_histogram(aes(x = value)) +
  facet_wrap(~ name, scales = "free") +
  labs(x = "", y = "Frequência") + 
  theme_classic()

ggsave("tex/imagens/indicadores.pdf", gg, width = 8, height = 2.75)
gg

IPTU %>% 
  filter(residencial == 1) %>% 
  ggplot(aes(x = verticalizacao, y = CA)) +
  geom_point(alpha = .75) +
  geom_abline(intercept = 0, slope = 1, color = "blue") +
  geom_smooth() +
  theme_classic() +
  coord_cartesian(xlim = c(0,50), ylim = c(0,20)) +
  theme(aspect.ratio = 1)

gg <- IPTU %>% 
  filter(residencial == 1, !is.na(ano_construcao), ano_construcao > 1950) %>% 
  mutate(ano_construcao = ano_construcao %>% round()) %>% 
  group_by(ano_construcao) %>% 
  # Cria quantis, o primeiro e último são o mínimo e máximo, respectivamente
  summarize(cota_parte = quantile(cota_parte,  0:6/6),
            CA = quantile(CA,  0:6/6),
            verticalizacao = quantile(verticalizacao,  0:6/6)) %>% 
  mutate(quantil = row_number() - 1,
         ano_construcao = as.Date(as.character(ano_construcao), format = "%Y")) %>% 
  pivot_longer(cota_parte:verticalizacao) %>% 
  group_by(ano_construcao, name) %>% 
  mutate(ymin = lag(value) %>% replace_na(0)) %>% 
  filter(quantil %in% 2:5) %>% 
  ggplot() +
  geom_ribbon(aes(x = ano_construcao, ymax = value, ymin = ymin, fill = factor(quantil)), color = "black") +
  facet_grid(rows = "name", scales = "free_y", 
             labeller = as_labeller(c("CA" = "CA", 
                                      "verticalizacao" = "Verticalização", 
                                      "cota_parte" = "Cota parte"))) +
  scale_x_date(limits = as.Date(c("1960", "2022"), format = "%Y")) +
  scale_fill_viridis_d(labels = c("20%", "40%", "60%", "80%")) +
  theme_bw() +
  labs(fill = "Quantil", x = "Ano da construção")

ggsave("tex/imagens/indicadores_tempo.pdf", gg, width = 8, height = 9)
ggsave("tex/imagens/indicadores_tempo_small.pdf", gg, width = 5.5, height = 6)

gg

Geosampa

GeoSampa é o portal de mapas e dados geoespaciais da cidade de São Paulo, mantido pela prefeitura. Ele fornece uma vasta gama de informações geográficas, incluindo mapas, dados demográficos, infraestruturas e muito mais. Este portal é uma ferramenta valiosa para pesquisadores, urbanistas e qualquer pessoa interessada em informações espaciais detalhadas sobre a cidade.

https://geosampa.prefeitura.sp.gov.br/

Base de Lotes: No GeoSampa, a base de lotes representa a divisão da cidade em pequenos segmentos, geralmente correspondentes a terrenos individuais ou condomínios. Esta base está organizada de forma que cada bairro tem seu próprio conjunto de dados e os dados de lotes para cada bairro podem ser baixados diretamente do site do GeoSampa em formato zip, contendo arquivos como .shp (shapefile), .dbf (database file), e .shx (index file).

# Extração do zip file com lotes de cada bairro
if (!"zip" %in% list.files(path = "dados/lotes")){
  for (file in list.files(path="dados/lotes/zip", full.names = FALSE) %>% 
       str_remove("\\.zip")){
    unzip(paste("dados/lotes/zip/", file, ".zip", sep = ""), 
          paste(file, "/", file, ".shp", sep = ""), 
          exdir = "dados/lotes/unzip")
    unzip(paste("dados/lotes/zip/", file, ".zip", sep = ""), 
          paste(file, "/", file, ".dbf", sep = ""), 
          exdir = "dados/lotes/unzip")
    unzip(paste("dados/lotes/zip/", file, ".zip", sep = ""), 
          paste(file, "/", file, ".shx", sep = ""), 
          exdir = "dados/lotes/unzip")
  }
}

# Junção dos shapes dos lotes de cada bairro em uma tabela
lotes <- list.files(path="dados/lotes/unzip", full.names = FALSE) %>% 
    paste("dados/lotes/unzip/", ., "/", ., ".shp", sep = "") %>% 
    lapply(read_sf) %>% 
    bind_rows %>% 
    st_set_crs("epsg:31983") 

Os lotes podem ser classificados de três formas

• Lotes fiscais: apresentam contribuintes do IPTU
• Espaço livre: espaços públicos
• Via de acesso: Ruas fechadas para acesso dentro de uma quadra

Os lotes fiscais podem ser unidades ou condomínios. Um prédio, por exemplo, fica em um único lote, mas dentro pode haver diversas unidades, então se configura como um condomínio. Não é possível saber através dos dados de lotes quantas unidades estão no condomínio, nem alguma forma de discriminá-los.

lotes %>% 
  mutate(condominio = case_when(lo_condomi == "00" ~ "Unidade",
                                lo_condomi != "00" ~ "Condomínio"),
         tipo = case_when(lo_tp_lote == "F" ~ "Fiscal",
                          lo_tp_lote == "M" ~ "Espaço livre",
                          lo_tp_lote == "V" ~ "Via de acesso",
                          .default = NA),
         area = st_area(geometry) %>% as.numeric()) %>% 
  st_drop_geometry() %>% 
  ggplot() +
  geom_violin(aes(x = factor(condominio), y = area, fill = tipo)) +
  scale_y_log10(labels = scales::comma_format(big.mark = ".")) +
  labs(title = "Área dos lotes em SP" , x = "", y = "Área em metros quadrados", fill = "Tipo de lote") +
  theme_classic()

Na base de dados de lotes, cada lote é identificado por três componentes principais:

1. Setor (lo_setor): Uma divisão maior dentro do município que agrupa várias quadras.
2. Quadra (lo_quadra): Uma subdivisão dentro de um setor que agrupa vários lotes.
3. Lote (lo_lote): A menor unidade de divisão, que representa um terreno ou uma parcela específica dentro de uma quadra. Essa estrutura de Setor-Quadra-Lote é crucial para identificar de forma única cada lote dentro do município. No código, esses identificadores são utilizados para manipular e cruzar os dados de lotes com outras bases de dados, como a base de IPTU.

lotes <- lotes %>% 
 filter(lo_tp_lote == "F") %>% # Seleção apenas de lotes fiscais
  mutate(lo_lote = ifelse(lo_lote == "0000", paste("CD", lo_condomi, sep = ""), lo_lote)) %>% 
  select(setor = lo_setor, quadra = lo_quadra, lote = lo_lote)

quadras <- read_sf("dados/quadras/SIRGAS_SHP_quadraMDSF.shp") %>% 
  st_set_crs("epsg:31983") %>% 
  filter(qd_tipo == "F") %>% # Apenas quadras fiscais
  select(setor = qd_setor, quadra = qd_fiscal) %>% 
  group_by(setor, quadra) %>% 
  summarize(geometry = st_union(geometry)) %>% 
  ungroup()

setores <- read_sf("dados/setor/SIRGAS_SHP_setorfiscal.shp") %>% 
  st_set_crs("epsg:31983") %>% 
  select(setor = st_codigo) %>% 
  group_by(setor) %>% 
  summarize(geometry = st_union(geometry)) %>% 
  ungroup()

Join IPTU - GeoSampa

# Join dos lotes com IPTU com base no SQL
IPTU.lote <- IPTU %>% 
  filter(residencial == 1) %>% 
  left_join(lotes, by = join_by(setor, quadra, lote)) %>% 
  ungroup()

IPTU.lote %>% 
  modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
unidades	545	0	2.4	16.4	1.0	1.0	2861.0
area_terreno	5736	0	248.8	1917.0	6.0	160.0	1442400.0
area_construida	12097	0	306.3	1954.9	2.0	120.0	400000.0
area_ocupada	3403	0	118.4	407.8	1.0	90.0	200000.0
pavimentos	51	0	1.8	1.8	1.0	2.0	52.0
ano_construcao	166	0	1981.8	14.6	1800.0	1981.0	2023.0
residencial	1	0	1.0	0.0	1.0	1.0	1.0
CA	100483	0	0.9	0.8	0.0	0.8	36.1
cota_parte	17359	0	208.3	1624.9	1.5	155.0	1442400.0
verticalizacao	53664	0	1.7	1.7	0.0	1.4	180.0

IPTU.quadra <- IPTU %>% 
  filter(residencial == 1) %>% 
  group_by(setor, quadra) %>% 
  summarize(unidades = sum(unidades),
            area_construida = sum(area_construida),
            verticalizacao = weighted.mean(verticalizacao, area_terreno),
            area_terreno = sum(area_terreno)) %>% 
  left_join(quadras, by = join_by(setor, quadra)) %>% 
  ungroup()

IPTU.quadra %>% 
  modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
unidades	976	0	77.6	160.6	1.0	34.0	4762.0
area_construida	15166	0	9976.7	19124.6	8.0	4813.0	521949.0
verticalizacao	33567	0	2.4	3.0	0.9	1.5	81.6
area_terreno	14846	0	8104.0	13914.7	50.0	6029.0	1447291.0

IPTU.setor <- IPTU %>% 
  filter(residencial == 1) %>% 
  group_by(setor) %>% 
  summarize(unidades = sum(unidades),
            area_construida = sum(area_construida),
            verticalizacao = weighted.mean(verticalizacao, area_terreno),
            area_terreno = sum(area_terreno)) %>% 
  left_join(setores, by = join_by(setor)) %>% 
  ungroup()

IPTU.setor %>% 
  modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
unidades	168	0	15827.8	9045.2	133.0	13694.0	46511.0
area_construida	168	0	2034480.2	1287662.5	16480.0	1682702.5	6822366.0
verticalizacao	168	0	3.5	2.2	1.1	2.7	12.5
area_terreno	168	0	1652592.8	1064452.3	30130.0	1410261.5	5087896.0

list(IPTU %>% filter(residencial == 1) %>% mutate(base = "bruta"), 
     IPTU.setor  %>% mutate(base = "Setor")  %>% filter(!geometry %>% st_is_empty()), 
     IPTU.quadra %>% mutate(base = "Quadra") %>% filter(!geometry %>% st_is_empty()), 
     IPTU.lote   %>% mutate(base = "Lote")   %>% filter(!geometry %>% st_is_empty())) %>% 
  lapply(function(x) x %>% 
           group_by(base) %>% 
           summarize(unidades = sum(unidades))) %>% 
  bind_rows %>% 
  pivot_wider(names_from = base, values_from = unidades) %>% 
  pivot_longer(2:4) %>% 
  mutate(missing = bruta - value,
         missing_percent = (100 *missing / bruta) %>% round(2) %>% paste("%")) %>% 
  select("Critério de join" = name,
         "Erros (unidades)" = missing,
         "Percentual de erro" = missing_percent) %>% 
  gt()

Critério de join	Erros (unidades)	Percentual de erro
Setor	0	0 %
Quadra	820	0.03 %
Lote	44319	1.67 %

gg <- IPTU.lote %>% 
  mutate(distancia_centro = st_distance(geometry, 
                                        read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% 
                                          st_set_crs("epsg:31983") %>% 
                                          filter(ds_nome == "SE")) %>% 
           as.numeric() %>% 
           cut(breaks = 10^3*(0:40), labels = FALSE)) %>% 
  st_drop_geometry() %>% 
  group_by(distancia_centro) %>% 
  summarize(verticalizacao = weighted.mean(verticalizacao, area_terreno)) %>% 
  ggplot(aes(x = distancia_centro, y = verticalizacao)) +
  geom_col() +
  geom_hline(yintercept = 1, linetype = "dotted") +
  labs(x = "Distância ao centro em quilômetros", 
       y = "Verticalização") + 
  scale_y_continuous(breaks = (1:6)) +
  theme_classic()

ggsave("tex/imagens/verticalizacao.pdf", gg, width = 6, height = 4)
gg

gg.lotes <- IPTU.quadra %>% 
  ggplot() +
  geom_sf(data = distrito, fill = "black") +
  geom_sf(aes(fill = verticalizacao, geometry = geometry), color = NA) +
  theme_void() +
  scale_fill_gradient(low = "#004D40FF", high = "#E0F2F1FF")

gg.setor <- IPTU.setor %>% 
  ggplot() +
  geom_sf(data = distrito, fill = "black") +
  geom_sf(aes(fill = verticalizacao, geometry = geometry), color = NA) +
  theme_void() +
  scale_fill_gradient(low = "#004D40FF", high = "#E0F2F1FF")

ggsave("tex/imagens/mapa_verticalizacao.pdf", 
       (gg.lotes | gg.setor), width = 20, height = 20)

(gg.lotes | gg.setor)

Join IPTU - Censo

O processo de cruzamento foi realizado com base na intersecção das geometrias dos setores censitários e dos lotes do IPTU. Cada setor censitário e cada lote do IPTU possui uma geometria associada, representando sua área geográfica no mapa. Ao cruzar essas geometrias, é possível identificar quais lotes estão contidos em cada setor censitário e vice-versa.

É importante destacar que, em casos onde um lote foi dividido entre dois ou mais setores censitários, ocorrerá uma intersecção em ambas as áreas. Para lidar com essa situação, foi calculado o percentual da área do lote que está contida em cada setor censitário.

# Join dados do IPTU com do Censo através da intersecção das geometrias
IPTU.censo <- censo %>% 
  st_intersection(IPTU.lote %>% st_as_sf()) %>% 
  as_tibble() %>% 
  rename(geometria_intersec = geometry) %>% 
  
  # Retomada das geometrias do setor e lotes
  left_join(censo %>% 
              as_tibble() %>% 
              select(id_setor_censitario, 
                     geometria_setor_censitario = geometry),
            by = join_by(id_setor_censitario)) %>% 
  left_join(IPTU.lote %>% 
              as_tibble() %>% 
              select(setor, quadra, lote, 
                     geometria_lote = geometry),
            by = join_by(setor, quadra, lote)) %>% 
  
  # Cálculo de quanto % do lote está dentro do setor
  mutate(percent_intersec = as.numeric(st_area(geometria_intersec) / st_area(geometria_lote))) 

IPTU.censo %>% 
  modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
v0001	1120	0	524.4	202.6	0.0	511.0	2221.0
v0002	543	0	227.8	87.1	0.0	222.0	1372.0
v0003	538	0	227.5	87.1	0.0	221.0	1371.0
v0004	24	0	0.2	1.0	0.0	0.0	170.0
v0005	11463	0	2.7	0.3	0.0	2.7	8.0
v0006	6421	0	13.4	10.3	0.0	11.2	100.0
v0007	471	0	194.7	72.7	0.0	190.0	950.0
area_setor	18531	0	52804.2	74529.5	196.3	41868.0	6945494.4
unidades	541	0	3.7	33.0	1.0	1.0	2861.0
area_terreno	5624	0	326.2	3606.1	6.0	162.0	1442400.0
area_construida	11985	0	460.5	4038.2	2.0	120.0	400000.0
area_ocupada	3350	0	140.3	1052.4	1.0	90.0	200000.0
pavimentos	51	0	1.9	2.2	1.0	2.0	52.0
ano_construcao	165	0	1981.7	14.6	1800.0	1981.0	2023.0
residencial	1	0	1.0	0.0	1.0	1.0	1.0
CA	99679	0	0.9	0.9	0.0	0.8	36.1
cota_parte	17099	0	240.4	2779.6	1.5	156.0	1442400.0
verticalizacao	53205	0	1.7	2.0	0.1	1.4	180.0
percent_intersec	285710	0	0.9	0.2	0.0	1.0	96.6

geometria <- IPTU.censo %>% 
                filter(id_setor_censitario == "355030810000143P") %>% 
                select(geometria_setor_censitario) %>% 
                st_as_sf()

zoom <- 60

bbox <- geometria %>% 
  st_transform("epsg:31983") %>% 
  st_bbox()

x_min <- (bbox$xmin - zoom) %>% as.numeric()
y_min <- (bbox$ymin - zoom) %>% as.numeric()
x_max <- (bbox$xmax + zoom) %>% as.numeric()
y_max <- (bbox$ymax + zoom) %>% as.numeric()

corte <- list(xmin = x_min - zoom, ymin = y_min - zoom, 
              xmax = x_max + zoom, ymax = y_max + zoom)

df <- IPTU.censo %>%
  st_as_sf() %>%
  st_transform("epsg:31983") %>% 
  st_crop(corte %>% unlist()) %>% 
  mutate(percent_intersec = round(percent_intersec * 100))


gg <- df %>% 
  mutate(texto = ifelse(percent_intersec != 100, percent_intersec, NA),
         color = ifelse(percent_intersec == 100, "blue", "red")) %>% 
  ggplot() +
  geom_sf(aes(geometry = geometria_setor_censitario), color = "black", fill = "grey", 
          lwd = 1.5, linetype ="solid") +
  geom_sf(aes(geometry = geometria_intersec, fill = color, color = color), 
          alpha = .25, linetype = "dashed", lwd = 1) +
  geom_sf_label(aes(geometry = geometria_intersec, label = texto), label.size = .01) +
  scale_fill_manual(values = c("red" = "red", "blue" = "blue")) +
  scale_color_manual(values = c("red" = "red", "blue" = "blue")) +
  labs(fill = "Decil de densidade") +
  theme_void() +
  theme(legend.position = "none") +
  coord_sf(xlim = c(x_min, x_max), ylim = c(y_min, y_max))
  

ggsave("tex/imagens/corte_lote.pdf", gg, width = abs(x_min - x_max)/40, height = abs(y_min- y_max)/40)

gg

Teoricamente o número de domicílios deveria ser próximo do número de unidades habitacionais enunciadas pelo IPTU. Caso haja mais unidades habitacionais do que domicílios, houve algum equívoco na medição do censo, visto que entre domicílios não contam apenas os ocupados. Caso haja mais domicílios, há unidades que não são contribuintes do IPTU. O que é preocupante é que casos em que o número é diferente não são exceções.

gg <- censo %>% 
  st_drop_geometry() %>% 
  select(id_setor_censitario, domicilios = v0002) %>% 
  left_join(IPTU.censo %>% 
              st_drop_geometry() %>% 
              filter(residencial == 1) %>% 
              group_by(id_setor_censitario) %>% 
              summarize(unidades = sum(unidades * percent_intersec))) %>% 
  mutate(unidades = replace_na(unidades, 0)) %>% 
  mutate(espectro_irregularidade = unidades / (unidades + domicilios)) %>% 
  ggplot() +
  geom_histogram(aes(x = espectro_irregularidade)) +
  labs(x = "Espectro de irregularidade: unidades / (unidades + domicilios)", y = "Frequência") +
  theme_classic() +
    theme(plot.caption = element_text(hjust = 0)) +
  scale_x_continuous(labels = scales::percent)


ggsave("tex/imagens/disparidade_censoIPTU.pdf", gg, width = 7, height = 3)
gg

Alguns setores censitários não encontram pares na base de lotes. Isso acontece principalmente por conta de loteamentos irregulares e favelas, que não são contribuintes do IPTU, e, portanto, não constam na base. Isso não prejudica a análise, pois estes loteamentos não dependem da regulamentação urbana, então mudanças nos instrumentos e indicadores não impactariam essas regiões.

Outras falhas decorrem do erro do join da base do IPTU com lotes, como apontado anteriormente, mas estes casos são negligenciáveis.

distrito <- read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% st_set_crs("epsg:31983")
lotes_irregulares <- read_sf("dados/lotes_irregulares/SIRGAS_SHP_loteamento.shp") %>%
  st_set_crs("epsg:31983")
favela <- read_sf("dados/favela/SIRGAS_SHP_favela.shp") %>% 
  st_set_crs("epsg:31983")


censo.pontos <- censo %>%
  select(id_setor_censitario, geometry, pessoas = v0001) %>%
  
  # Um ponto a cada 1.000 pessoas
  mutate(pontos = pessoas %/% 1000 + (runif(length(pessoas)) * 1000 < pessoas %% 1000)) %>% 
  select(-pessoas) %>% 
  as_tibble() %>% 
  left_join(
            # Seleção dos setores censitários que não apresentam nenhum contribuinte do IPTU
            censo %>% 
              anti_join(IPTU.censo) %>% 
              st_drop_geometry() %>% 
              select(id_setor_censitario) %>% 
              mutate(erro = TRUE)) %>% 
  mutate(erro = replace_na(erro, FALSE)) %>% 
  filter(pontos > 0) %>% 
  mutate(samples = purrr::map2(geometry, pontos, ~st_sample(.x, size = .y))) %>%
  unnest(cols = samples) %>%
  as_tibble()

gg <- censo.pontos %>% 
  ggplot() +
  geom_sf(data = distrito, color = NA) +
  geom_sf(data = st_union(favela %>% st_union() %>% st_buffer(10),
                          lotes_irregulares %>% st_union() %>% st_buffer(10)),
          aes(fill = "Favelas e lotes irregulares"), 
        color = "black", size = .1, alpha = .5) +
  geom_sf(aes(geometry = samples, color = erro), alpha = .5, size = .8) +
  scale_color_manual(values = c("FALSE" = NA,#248232
                               "TRUE" = "#D32934FF"),
                     labels = NULL) +
  scale_fill_manual("", values = c("Favelas e lotes irregulares" = "#2BAA92FF")) +
  labs(title = "<span style='font-size: 34pt;'>População em <span style = 'color:#D32934FF;'>áreas sem registro de IPTU</span> geralmente estão em <br><span style = 'color:#2BAA92FF;'>favelas ou lotes irregulares</span> (cada ponto representa 1000 pessoas)</span>") + 
  theme_void() +
  theme(plot.title = element_markdown(), legend.position = "none")

# scale_colour_paletteer_d("lisa::AndyWarhol_2")

ggsave("tex/imagens/mapa_pontos.pdf", gg, width = 16, height = 20)
gg

Análise de resultados via join

df <- IPTU.censo %>% 
  filter(residencial == 1) %>% 
  
  st_drop_geometry() %>% 
  
  group_by(id_setor_censitario) %>% 
  summarize(# Setor censitário
            populacao = median(v0001),
            area_setor = median(area_setor) %>% as.numeric(),
            domicilios = median(v0002),
            domicilios_ocupados = median(v0007),
            
            # Lote
            unidades = sum(unidades * percent_intersec),
            area_terreno = sum(area_terreno * percent_intersec),
            area_construida = sum(area_construida * percent_intersec),
            area_ocupada = sum(area_ocupada * percent_intersec),
            ) %>% 
  
  mutate(densidade = populacao * 10 ^ 6 / area_setor,
         cota_parte = area_terreno / unidades,
         CA = area_construida / area_terreno,
         verticalizacao = area_construida / area_ocupada,
         espectro_irregularidade = unidades / (unidades + domicilios))

modelsummary::datasummary_skim(df)

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
populacao	1120	0	420.0	213.4	0.0	397.0	2221.0
area_setor	18531	0	37267.5	109998.2	196.3	24549.0	6945494.4
domicilios	543	0	189.6	88.1	0.0	182.0	1372.0
domicilios_ocupados	471	0	162.5	75.3	0.0	156.0	950.0
unidades	18524	0	141.1	97.1	0.0	130.1	1343.0
area_terreno	18531	0	14702.9	18405.7	0.0	11595.9	1445351.6
area_construida	18531	0	18094.3	13503.4	0.0	15739.1	163349.7
area_ocupada	18531	0	7008.5	6150.2	0.0	5626.0	64738.0
densidade	18341	0	28519.9	36697.1	0.0	17034.3	1162304.8
cota_parte	18156	0	1282.3	14374.8	1.8	133.4	516201.6
CA	18238	0	2.3	2.4	0.0	1.0	27.5
verticalizacao	17941	0	4.7	5.5	1.0	1.8	77.4
espectro_irregularidade	18360	0	0.4	0.2	0.0	0.4	1.0

df %>% 
  select(id_setor_censitario, populacao, domicilios, unidades, area_setor) %>% 
  bind_rows(censo %>% 
              st_drop_geometry() %>% 
              select(id_setor_censitario, populacao = v0001, domicilios = v0002, area_setor) %>% 
              mutate(unidades = 0) %>% 
              anti_join(df %>% select(id_setor_censitario))) %>% 
  filter(populacao > 0) %>% 
  mutate(densidade = populacao/area_setor,
         taxa_irregular = (domicilios - unidades) / domicilios,
         espectro_irregularidade = unidades / (unidades + domicilios)) %>% 
  lm(log(densidade) ~ espectro_irregularidade, .) %>% 
  summary()

## 
## Call:
## lm(formula = log(densidade) ~ espectro_irregularidade, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.0873  -0.4522   0.0812   0.7049   4.1956 
## 
## Coefficients:
##                         Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)             -4.10609    0.01168 -351.450  < 2e-16 ***
## espectro_irregularidade  0.14028    0.03279    4.278 1.89e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.224 on 26833 degrees of freedom
## Multiple R-squared:  0.0006816,  Adjusted R-squared:  0.0006443 
## F-statistic:  18.3 on 1 and 26833 DF,  p-value: 1.893e-05

modelo1 <- df %>% 
  lm(densidade ~ cota_parte + verticalizacao + CA, data = .) %>% 
  summary()

modelo2 <- df %>% 
  filter(abs(espectro_irregularidade - .5) < .25) %>% 
  lm(densidade ~ cota_parte + verticalizacao + CA, data = .) %>% 
  summary()

modelo3 <- df %>% 
  filter(abs(espectro_irregularidade - .5) < .1) %>% 
  lm(densidade ~ cota_parte + verticalizacao + CA, data = .) %>% 
  summary()

modelsummary::modelsummary(
  list("Todos setores" = modelo1, "Balanço entre 25% e 75%" = modelo2, "Balanço entre 40% e 60%" = modelo3), 
  gof_omit = "RMSE", estimate = "{estimate} {stars}", 
  statistic = "({std.error})")

	Todos setores	Balanço entre 25% e 75%	Balanço entre 40% e 60%
(Intercept)	8603.264 ***	4242.744 ***	348.108
	(311.878)	(347.853)	(608.823)
cota_parte	0.095 ***	-0.018	-4.176 *
	(0.016)	(0.057)	(1.712)
verticalizacao	1058.806 ***	1146.076 ***	1283.914 ***
	(59.059)	(61.902)	(76.170)
CA	6524.131 ***	7221.037 ***	7874.627 ***
	(132.933)	(138.888)	(171.647)
Num.Obs.	18531	15616	9363
R2	0.314	0.374	0.410
R2 Adj.	0.314	0.374	0.409

split <- rsample::initial_split(df, prop = 0.7)
train <- rsample::training(split)
test <- rsample::testing(split)

rf_model <- ranger::ranger(densidade ~ ., 
                     data = train, num.trees = 5000, importance = "permutation")

## Growing trees.. Progress: 100%. Estimated remaining time: 0 seconds.

var.importance <- ranger::importance(rf_model)
gg <- ggplot(NULL, aes(x = var.importance,
                 y = fct_reorder(names(var.importance) %>% str_replace_all("_", " "),
                                 var.importance))) +
    geom_col(fill = "darkblue", alpha = 0.75) +
    labs(x = "Importância", y = "", title = "Variáveis mais relevantes para explicar a densidade", 
         subtitle = "Resultados da metodologia via join") +
    theme_classic()

ggsave("tex/imagens/importancia.pdf", gg, width = 7, height = 4)

remove(rf_model)

gg

Raster

geoms.IPTU <- IPTU.lote %>% 
  filter(! st_is_empty(geometry)) %>% 
  st_as_sf() %>% 
  mutate(area = st_area(geometry) %>% as.numeric()) %>% 
  select(id = setor, area, unidades, area_construida, area_ocupada, area_terreno)

geoms.censo <- censo %>% 
  mutate(area = st_area(geometry) %>% as.numeric()) %>% 
  select(id = id_setor_censitario, area, populacao = v0001, domicilios = v0002)

bbox <- st_bbox(geoms.censo)

raster <- raster::raster(xmn = bbox["xmin"], xmx = bbox["xmax"], ymn = bbox["ymin"], ymx = bbox["ymax"],
            res = 800, crs = st_crs(geoms.censo)$proj4string, vals = NA)

result.IPTU <- exact_extract(raster, geoms.IPTU, include_xy = T, 
                             include_cols = c("id", "area", "unidades", "area_construida", "area_ocupada",
                                              "area_terreno"), 
                             force_df = T, coverage_area = T, progress = F) %>% 
  bind_rows %>% 
  mutate(percent_intersect = coverage_area / area,
         across(unidades:area_terreno, ~ .x * percent_intersect)) %>% 
  group_by(x, y) %>% 
  summarize(unidades = sum(unidades),
            area_ocupada = sum(area_ocupada),
            area_construida = sum(area_construida),
            area_terreno = sum(area_terreno)) %>% 
  ungroup() %>% 
  mutate(cota_parte = area_terreno / unidades,
         CA = area_construida / area_terreno,
         verticalizacao = area_construida / area_ocupada)

result.censo <- exact_extract(raster, geoms.censo, include_xy = T, 
                        include_cols = c("id", "area", "populacao", "domicilios"), 
                        force_df = T, coverage_area = T, progress = F) %>% 
  bind_rows %>% 
  mutate(percent_intersect = coverage_area / area,
         across(c(area, populacao, domicilios), ~ .x * percent_intersect)) %>% 
  group_by(x, y) %>% 
  summarize(populacao = sum(populacao),
            area = sum(area),
            domicilios = sum(domicilios)) %>% 
  mutate(densidade = populacao * 10^6 / area)

result <- full_join(result.IPTU, result.censo) %>% 
  mutate(espectro_irregularidade = unidades / (unidades + domicilios))

result %>% 
  modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
x	59	0	331929.5	11037.7	313794.9	329794.9	360194.9
y	91	0	7383596.5	18725.5	7343788.7	7388588.7	7415788.7
unidades	1256	52	2083.5	2099.4	0.0	1733.9	21056.9
area_ocupada	1256	52	103486.6	64038.7	1.5	113195.8	251213.6
area_construida	1256	52	267192.3	267862.6	1.5	216295.9	2321166.4
area_terreno	1256	52	217116.1	125745.5	7.8	238577.3	840996.4
cota_parte	1254	52	2647.1	42935.6	6.5	140.7	1415766.7
CA	1255	52	1.3	1.2	0.0	0.9	9.6
verticalizacao	1246	52	2.6	2.3	1.0	1.8	28.9
populacao	2403	0	4367.6	4804.2	0.0	2620.9	29598.3
area	1787	0	580109.1	157947.2	0.1	640000.0	640000.0
domicilios	2407	0	1905.6	2154.4	0.0	1108.2	17874.7
densidade	2355	0	7338.6	7859.5	0.0	4962.7	46247.3
espectro_irregularidade	1256	52	0.4	0.2	0.0	0.4	1.0

# r <- raster::rasterize(x = result %>% select(x, y), y = raster, 
#                        field = result %>% 
#                          select(CA, cota_parte, verticalizacao, densidade) %>% 
#                          mutate(across(c(CA, cota_parte, verticalizacao), ~ replace_na(.x, 0))))

gg <- result %>% 
  pivot_longer(c(CA, cota_parte, verticalizacao, densidade)) %>% 
  group_by(name) %>% 
  mutate(value = ntile(value, 10) %>% as.factor(),
         name = case_when(name == "CA" ~ "Coeficiente de Aproveitamento (CA)",
                          name == "cota_parte" ~ "Cota-parte",
                          name == "densidade" ~ "Densidade populacional",
                          name == "verticalizacao" ~ "Verticalização",
                          .default = NA)) %>% 
  ggplot() +
  geom_sf(data = distrito) +
  geom_tile(aes(x = x, y = y, fill = value), alpha = .75, color = "grey") +
  # facet_wrap(~name, ncol = 4) +
  theme_void() +
  theme(strip.text = element_text(size = 12)) +
  labs(fill = "Decil") +
  scale_fill_viridis_d()

ggsave("tex/imagens/rasters.pdf", gg + facet_wrap(~name, ncol = 2), width = 8, height = 12)
ggsave("tex/imagens/rasters_wide.pdf", gg + facet_wrap(~name, ncol = 4), width = 16, height = 6)

gg + facet_wrap(~name, ncol = 2)

(result %>%
  filter(area > max(area) * .99) %>%
  arrange(desc(densidade)) %>%
  mutate(top = ifelse(row_number() <= 5, row_number(), NA) %>% as.factor()) %>% 
  ggplot() +
  geom_sf(data = distrito, color = "grey", aes(text = ds_nome)) +
  geom_tile(aes(x = x, y = y, fill = top, text = paste("pop: ", populacao, "\n", 
                                                       "unidades: ", unidades, "\n",
                                                       "domicilios: ", domicilios, "\n",
                                                       "densidade: ", densidade, "\n",
                                                       "area: ", area,
                                                       sep = ""))) +
  scale_fill_viridis_d() +
  theme_void()) %>%
  plotly::ggplotly()

tabela <- result %>%
  filter(area > max(area) * .99) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>% 
  arrange(desc(densidade)) %>%
  head(5) %>% 
  mutate(localizacao = c("1. Paraisópolis", "2. Heliópolis ", "3. Sé (Bela Vista)", "4. Paraisópolis", "5. Heliópolis")) %>% 
  select("Localização" = localizacao, "População" = populacao, "Domicílios (Censo)" = domicilios,
         "Unidades (IPTU)" = unidades, "Espectro irregularidade" = espectro_irregularidade,
         "Densidade habitacional" = densidade, "Área" = area) %>% 
  pivot_longer(cols = -Localização, names_to = "Variável") %>% 
  pivot_wider(id_cols = Variável, names_from = Localização, values_from = value) %>% 
  gt() %>% 
  tab_spanner("Favelas", columns = c(2,3,5,6)) %>% 
  fmt_number(rows = c(1,2,3,5,6), decimals = 0, sep_mark = ".", dec_mark = ",") %>% 
  fmt_percent(rows = 4) 

gtsave(tabela %>% tab_options(table.font.size = 12), "tex/tabelas/top5.tex")
gtsave(tabela %>% tab_options(table.font.size = 9), "tex/tabelas/top5_small.tex")

tabela

Variável	Favelas				3. Sé (Bela Vista)
	1. Paraisópolis	2. Heliópolis	4. Paraisópolis	5. Heliópolis
População	29.598	25.280	23.824	22.920	24.576
Domicílios (Censo)	11.655	10.178	9.361	9.001	17.875
Unidades (IPTU)	0	1.857	7	3	21.057
Espectro irregularidade	0.00%	15.43%	0.08%	0.03%	54.09%
Densidade habitacional	46.247	39.500	37.225	35.813	38.400
Área	640.000	640.000	640.000	640.000	640.000

gg <- result %>% 
  drop_na() %>% 
  ggplot() +
  geom_sf(data = distrito, color = NA) +
  geom_tile(aes(x = x, y = y, fill = espectro_irregularidade), alpha = .75, color = "grey") +
  theme_void() +
  scale_fill_gradient2(high = "#2F191B", low = "#D32934", mid = "white", na.value = NA,
                       midpoint = .5, limits = c(0,1)) +
  theme(legend.position = c(.8, .35), legend.title = element_markdown(size = 25), 
        plot.title = element_markdown(), legend.text = element_text(size = 25),
        legend.key.size = unit(1.75, 'cm')) +
  labs(fill = "Espectro de irregularidade", title = "<span style='font-size: 35pt;'>Áreas com <span style = 'color:#D32934;'>menos domicílios registrados no IPTU, comparados <br>ao censo</span> geralmente estão em regiões afastadas do centro</span>")

ggsave("tex/imagens/balanco_raster.pdf", gg, width = 16, height = 20)

gg

Análise resultados via Raster

modelo1 <- result %>% 
  lm(densidade ~ cota_parte + verticalizacao + CA, data = .) %>% 
  summary()

modelo2 <- result %>% 
  filter(populacao > 0) %>% 
  lm(log(densidade) ~ cota_parte + verticalizacao + CA, data = .) %>% 
  summary()

modelo3 <- result %>% 
  filter(abs(espectro_irregularidade - .5) < .1) %>% 
  lm(densidade ~ cota_parte + verticalizacao + CA, data = .) %>% 
  summary()

modelo4 <- result %>% 
  filter(abs(espectro_irregularidade - .5) < .1,
         populacao > 0) %>% 
  lm(log(densidade) ~ cota_parte + verticalizacao + CA, data = .) %>% 
  summary()

modelos <- list("Nível (A)" = modelo1,
                "Log (B)" = modelo2,
                "Nível (C)" = modelo3,
                "Log (D)" = modelo4)
  
tabela <- modelsummary(modelos,
             estimate = c("Coeficiente" = "estimate"),
             statistic = c("*" = "stars"),
             shape = term ~ model + statistic,
             output = "gt",
             align = "lrlrlrlrl",
             add_rows = modelos %>% 
               sapply(get_gof) %>% 
               rbind %>% 
               as_tibble(rownames = 'variable') %>% 
               unnest(cols = c(`Nível (A)`, `Log (B)`, `Nível (C)`, `Log (D)`)) %>% 
               mutate(across(c(`Nível (A)`, `Log (B)`, `Nível (C)`, `Log (D)`), ~ 
                               ifelse(variable == "nobs", round(.) %>% as.character(), sprintf('%.3f',.))),
                      variable = case_when(variable == "r.squared" ~ "R2",
                                           variable == "adj.r.squared" ~ "R2 ajustado",
                                           variable == "nobs" ~ "Observações",
                                           .default = NA)) %>% 
               drop_na() %>% 
               mutate(!!!list("As" = "", "Bs" = "", "Cs" = "", "Ds" = "")) %>% 
               select(variable, "Nível (A)", As, "Log (B)", Bs, "Nível (C)", Cs, "Log (D)", Ds,)) %>% 
  tab_spanner("Todos os setores", columns = 2:5) %>% 
  tab_spanner("Espectro irregularidade: 40 a 60%", columns = 6:9) %>% 
  tab_style(style = cell_borders(sides = "top",  weight = px(0.5)),
            locations = cells_body(rows = 5)) 

gtsave(tabela %>% tab_options(table.font.size = 14), "tex/tabelas/regressao.tex")
gtsave(tabela %>% tab_options(table.font.size = 12), "tex/tabelas/regressao_small.tex")

tabela

	Todos os setores				Espectro irregularidade: 40 a 60%
	Nível (A)		Log (B)		Nível (C)		Log (D)
	Coeficiente	*	Coeficiente	*	Coeficiente	*	Coeficiente	*
(Intercept)	11046.375	***	8.987	***	8632.310	***	9.109	***
cota_parte	-0.004		0.000	***	-6.973	***	-0.002	***
verticalizacao	-249.059		-0.075	**	-48.242		-0.026
CA	1211.554	***	0.243	***	1672.931	***	0.143	***
R2	0.022		0.052		0.249		0.328
R2 ajustado	0.020		0.050		0.245		0.325
Observações	1255		1254		553		553

split <- rsample::initial_split(result %>% 
                                  drop_na() %>% 
                                  filter(abs(espectro_irregularidade - .5) <.1,
                                         area > 63950) %>% 
                                  select(-c(populacao, domicilios, espectro_irregularidade, area)), 
                                prop = 0.7)
train <- rsample::training(split)
test <- rsample::testing(split)

rf_model <- ranger::ranger(densidade ~ ., 
                     data = train, num.trees = 5000, importance = "permutation")
rf_preds <- predict(rf_model, data = test)$predictions

rf_model_restrito <- ranger::ranger(densidade ~ CA + cota_parte + verticalizacao, 
                     data = train, num.trees = 5000, importance = "permutation")
rf_preds_restrito <- predict(rf_model_restrito, data = test)$predictions

gg <- ggplot(NULL, aes(x = rf_model$variable.importance,
                 y = fct_reorder(names(rf_model$variable.importance),
                                 rf_model$variable.importance))) +
    geom_col(fill = "darkblue", alpha = 0.75) +
    labs(x = "Importância", y = "Variável", title = "") +
    theme_classic()

ggsave("tex/imagens/var_importance.pdf", gg, width = 6.5, heigh = 4)
gg

rtree <- rpart::rpart(densidade ~ CA + cota_parte + verticalizacao, data = train)

visNetwork::visTree(rtree, direction = "LR", collapse = TRUE, legend = FALSE)

lm_model <- lm(densidade ~ CA + cota_parte + verticalizacao, data = train)
lm_preds <- predict(lm_model, newdata = test)

1 - sum((rf_preds - test$densidade)^2)/sum((test$densidade - mean(test$densidade))^2)

## [1] 0.7979854

1 - sum((rf_preds_restrito - test$densidade)^2)/sum((test$densidade - mean(test$densidade))^2)

## [1] 0.4571126

1 - sum((lm_preds - test$densidade)^2)/sum((test$densidade - mean(test$densidade))^2)

## [1] 0.1076516

split <- rsample::initial_split(result %>% 
                                  drop_na() %>% 
                                  mutate(densidade = populacao / area_terreno) %>% 
                                  filter(abs(espectro_irregularidade - .5) <.1,
                                         area > 63950) %>% 
                                  select(-c(populacao, domicilios, espectro_irregularidade, area)), 
                                prop = 0.7)
train <- rsample::training(split)
test <- rsample::testing(split)

rf_model <- ranger::ranger(densidade ~ ., 
                     data = train, num.trees = 5000, importance = "permutation")
rf_preds <- predict(rf_model, data = test)$predictions

gg <- ggplot(NULL, aes(x = rf_model$variable.importance,
                 y = fct_reorder(names(rf_model$variable.importance),
                                 rf_model$variable.importance))) +
    geom_col(fill = "darkblue", alpha = 0.75) +
    labs(x = "Importância", y = "Variável", title = "") +
    theme_classic()

ggsave("tex/imagens/var_importance_densmod.pdf", gg, width = 6.5, heigh = 4)
gg